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ABSTRACT 

Computational  micromagnetics  in  three  dimensions  is  of  increasing  interest  with  the  development  of  magnetostrictive 
sensors  and  actuators.  In  solving  the  Landau-Lifshitz-Gilbert  (LLG)  equation,  the  governing  equation  of  magnetic 
dynamics  for  ferromagnetic  materials,  we  need  to  evaluate  the  effective  field.  The  effective  field  consists  of  several 
terms,  among  which  the  demagnetizing  field  is  of  long-range  nature.  Evaluating  the  demagnetizing  field  directly 
requires  work  of  0(N 2)  for  a  grid  of  N  cells  and  thus  it  is  the  bottleneck  in  computational  micromagnetics.  A 
fast  hierarchical  algorithm  using  multipole  approximation  is  developed  to  evaluate  the  demagnetizing  field.  We 
first  construct  a  mesh  hierarchy  and  divide  the  grid  into  boxes  of  different  levels.  The  lowest  level  box  is  the 
whole  grid  while  the  highest  level  boxes  are  just  cells.  The  approximate  field  contribution  from  the  cells  contained 
in  a  box  is  characterized  by  the  box  attributes,  which  are  obtained  via  multipole  approximation.  The  algorithm 
computes  field  contributions  from  remote  cells  using  attributes  of  appropriate  boxes  containing  those  cells,  and  it 
computes  contributions  from  adjacent  cells  directly.  Numerical  results  have  shown  that  the  algorithm  requires  work 
of  O(NlogN)  and  at  the  same  time  it  achieves  high  accuracy.  It  makes  micronragnetic  simulation  in  three  dimensions 
feasible. 

Keywords:  Micromagnetics,  demagnetizing  field,  multipole  approximation 

1.  INTRODUCTION 

Computational  micromagnetics  in  three  dimensions  is  of  increasing  interest  with  the  development  of  magnetostrictive 
sensors  and  actuators.  The  Landau-Lifshitz-Gilbert  (LLG)  equation  is  the  governing  equation  of  magnetic  dynamics 
for  ferromagnetic  materials.1  By  integrating  the  LLG  equation,  we  can  solve  for  the  evolution  of  magnetization 
as  well  as  the  steady  state  magnetization  profile.  This  helps  us  understand  the  underlying  physical  principles, 
characterize  material  properties,  and  design  and  control  magnetostrictive  transducers. 

The  effective  magnetic  field  Heff  needs  to  be  evaluated  in  solving  the  LLG  equation.  Heff  is  the  sum  of  several 
terms:  the  externally  applied  field  Hext,  the  anisotropy  field  Hanis  due  to  the  crystalline  anisotropy,  the  exchange 
field  Hexch  due  to  the  quantum-mechanical  exchange  effect  between  nearest  neighbours,  and  the  demagnetizing 
field  Hdemag  produced  by  the  whole  magnetization  distribution.  The  last  term  is  non-local,  because  the  decay  of 
the  demagnetizing  field  with  distance  is  so  slow  that  all  interactions  must  be  accounted  for.  Therefore  evaluation 
of  Hdemag  is  the  most  time-consuming  part  and  thus  the  bottleneck  in  computational  micromagnetics.  For  a 
discretization  grid  of  N  cells,  an  amount  of  work  of  0(N2)  is  required  to  evaluate  all  pairwise  interactions.  In  three 
dimensional  micromagnetics,  to  be  of  physical  interest,  the  simulation  is  usually  involved  in  thousands  of  cells.  As  a 
result,  the  computation  would  be  prohibitive  if  we  calculate  the  demagnetizing  field  directly. 

Many  techniques  have  been  proposed  to  speed  up  the  evaluation  of  Hdemag-  The  simplest  method  is  truncation 
of  the  interaction  range.2  Although  it  reduces  the  computation  time  to  O(N),  the  loss  of  accuracy  is  significant. 
The  Fast  Fourier  Transform  (FFT)  technique  is  used  widely  and  it  requires  work  of  0(NlogN ).3~5  Multipole 
approximation  is  another  useful  method  to  accelerate  the  computation  of  Hdemag- 6  4  The  general  strategy  of  the 
multipole  algorithm  is  clustering  cells  at  different  spatial  scales  and  using  multipole  expansions  to  evaluate  the 

*  This  research  was  supported  by  the  Army  Research  Office  under  the  ODDR&E  MURI97  Program  Grant  No.  DAAG55-97-1-0114 
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interactions  between  clusters  that  are  sufficiently  far  away  from  each  other.  The  interactions  between  nearby  cells 
are  calculated  directly.  Greengard  adopted  a  two-pass  (upward  and  downward)  nrultipole  algorithm  in  electrostatic 
calculations  and  obtained  O(N)  complexity.7  Following  this  approach,  a  fast  algorithm  fully  implementing  the 
multipole  and  local  expansions  of  the  field  integral  was  shown  to  yield  O(N)  computation  time  in  2D  nricromagnetics.4 
But  this  technique  can’t  be  generalized  to  3D  directly  because  the  closed  forms  of  multipole  and  local  expansions  in 
3D  do  not  exist.  Blue  and  Scheinfein  used  only  a  single  upward  pass  in  multipole  expansion  for  2D  micromagnetics, 
which  yielded  a  computation  time  of  O(NlogN).6  Direct  generalization  of  this  algorithm  to  3D  case  is  not  trivial 
either. 

Inspired  by  the  work  of  Blue  and  Scheinfein,  this  paper  is  aimed  at  proposing  a  simple  but  efficient  hierarchical 
algorithm  for  evaluation  of  Hdemag  in  3D.  We  first  construct  a  mesh  hierarchy  and  divide  the  grid  into  boxes  of 
different  levels.  The  approximate  field  contribution  from  the  cells  contained  in  a  box  is  characterized  by  the  box 
attributes,  which  are  obtained  via  multipole  approximation.  The  algorithm  computes  field  contributions  from  remote 
cells  using  attributes  of  appropriate  boxes  containing  those  cells,  and  it  computes  contributions  from  adjacent  cells 
directly.  The  paper  is  organized  as  follows:  Section  2  describes  the  fast  algorithm  in  detail,  Section  3  presents  the 
numerical  results  and  Section  4  provides  the  conclusions. 


2.  FAST  HIERARCHICAL  ALGORITHM  USING  MULTIPOLE  APPROXIMATION 


Assuming  that  a  ferromagnetic  body  with  magnetization  M  occupies  a  region  V,  the  demagnetizing  field  at  r, 
Hdemag(r),  can  be  written  down  directly  from  magnetostatics8: 
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r')V  •  M(r')  ,  , 

- rr - dr 

r-r'  1 1 3 
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r')M(r')  •  n(r') 
||  r-r'  1 1 3 
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(1) 


where  dr,,ds/  are  the  volume  element  and  the  area  element,  respectively,  the  integrals  are  taken  over  the  space 
variable  r',  “V-”  is  the  divergence  operator  with  respect  to  r',  j|  •  ||  is  the  Euclidean  norm  of  a  3D  vector,  dV  is  the 
boundary  of  V,  and  n(r')  is  the  unit  outward  normal  at  r'.  It  is  clear  from  (1)  that  Hdemag  depends  on  the  whole 
distribution  of  M. 


In  this  paper,  we  consider  a  ferromagnetic  body  with  rectangular  geometry.  We  discretize  the  body  into  cells 
of  size  do-  The  numbers  of  cells  along  the  x,y,  z  axes  are  denoted  by  Nx.  Ny  and  Nz,  so  the  total  number  of  cells 
N  =  NxNyNz.  M  is  assumed  to  be  constant  within  each  cell  (including  at  the  boundary). 

In  this  section,  we  will  first  look  at  the  demagnetizing  field  Hd  produced  by  a  cell  with  magnetization  M.  The 
analytic  expression  for  Hd  turns  out  to  be  very  complicated.  We  then  show  that  Hd  can  be  approximated  by  the 
field  produced  by  a  magnetic  dipole  with  moment  Mc?q  located  at  the  center  of  the  cell.  Although  the  dipole  formula 
looks  much  simpler,  direct  pairwise  evaluation  is  still  very  time-consuming  when  N  is  large.  A  hierarchical  algorithm 
incorporating  nrultipole  approximation  is  then  used  to  accelerate  the  evaluation  of  Hdemag- 


2.1.  Demagnetizing  Field  Contribution  from  a  Cell 

We  want  to  evaluate  the  demagnetizing  field  Hd  at  the  center  of  cell  ( i,j,k )  produced  by  cell  (: i',j',k ')  with 
magnetization  M.  Denote  the  region  that  cell  ( ')  occupies  as  P.  Since  M  is  constant  within  P,  we  have 
V  •  M  =  0  and  only  the  second  term  in  (1)  survives.  Let  Hdx,  Hdy,  Hdz  be  the  x,y,z  components  of  Hd,  and  let 
nx  =  i  —  i' ,  ny  =  j  —  j' ,  nz  =  k  —  k! .  It  follows  that 


where  (x  y  z)T  =  do(nx 
expression  for  Hdx, 


didx  — 
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dx'dy'dz' , 
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ion  [y/ (x  -  x')2  +  (y  —  y')2  +  [z  -  z')2]3 
ny  nz)T ,  and  x',  y',  z'  G  [—  4r,  4j*].  Carrying  out  the  integral,  we  will  get  the  analytic 
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where  Mx,  My ,  Mz  are  the  x,  y ,  z  components  of  M,  and  lXi,  lyi,  lzz,  i  =  1,  ■  •  •  8  are  functions  of  nx,  ny,  nz,  as  defined 
in  Appendix  A.  Similar  expressions  can  be  obtained  for  Hdy  and  Hdz-  Note  that  Hd  does  not  depend  on  the  cell 
size  d0;  instead,  it  depends  only  on  M  and  ( nx  nv  nz)T . 


2.2.  Approximation  to  Hj  by  the  Dipole  Formula 


Expanding 


in  (2),  we  will  get 
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where  r=  ( x  y  z)T .  Plugging  the  first  order  approximation  in  (5)  rather  than  (4)  into  (2)  and  carrying  out  the 
integral  leads  to  the  approximation  to  Hdx,  denoted  as  Hdx. 


Hdx  — 


[3(M  •  r)x  -  Mx  ||  r  ||2]dg 


(6) 


Similarly,  we  can  get  Hdv  and  Hciz.  It’s  easy  to  see  that,  letting  Hd  =  (Hdx  Hdy  Hdz)T  and  r  =  r/  ||  r  ||, 


*  [3(M  •  ?)r  -  M ]dg 

Hd - FF - 

which  is  the  familiar  dipole  formula. 

Defining  the  error 

||  Hd  -  Hd  || 

6  II  Hd  ||  ’ 


(7) 

(8) 


we  expect  e  to  be  of  order  0(l/(n2  +  n2  +  n2))  from  (2),  (5)  and  x' ,  y',  z’  €  [—  4^,  4p]. 

The  dipole  formula  (7)  is  much  simpler  than  (3),  which  we  shall  refer  as  the  integral  formula.  But  we  can’t  be 
too  optimistic  here.  When  N  goes  large,  the  computation  time  of  direct  pairwise  evaluation  increases  with  0(N2). 
In  other  words,  halving  the  cell  size  will  multiply  the  computation  time  by  64.  Even  with  the  dipole  formula,  the 
computation  may  become  prohibitive  well  before  the  grid  is  fine  enough.  Thus  we  resort  to  the  following  hierarchical 
algorithm. 


2.3.  A  Hierarchical  Evaluation  Algorithm 

As  pointed  out  earlier,  the  cell  size  do  does  not  appear  in  the  analytic  expression  of  Hd,  therefore  from  now  on  we 
will  assume  unit  cell  size.  We  start  with  a  3D  grid  of  N  cells  and  embed  the  grid  into  a  box  of  size  2m,  with  the 
smallest  possible  m  to  enclose  the  grid.  Mesh  level  0  corresponds  to  the  entire  computational  box,  while  mesh  level 
L  +  1  is  obtained  from  level  L  by  subdividing  each  box  into  eight  subboxes.  Continue  this  process  until  there  is  at 
most  one  cell  in  each  box.  A  tree  structure  is  imposed  on  this  mesh  hierarchy,  so  that  if  B  is  a  box  at  level  L,  the 
eight  boxes  at  level  L  +  1  obtained  by  subdividing  B  are  considered  its  children. 

As  we  will  show  in  the  next  subsection,  the  approximate  demagnetizing  field  contribution  from  cells  contained 
in  each  box  can  be  characterized  with  two  attributes  of  that  box.  Assume  that  we  have  got  these  attributes  for  all 
boxes  and  pick  a  threshold  value  (3.  The  following  steps  can  be  taken  to  evaluate  the  field  at  the  center  of  each  cell 
in  the  grid: 


•  Step  1.  Put  all  boxes  of  level  1  on  a  stack  (since  the  box  of  level  0  will  contain  the  current  cell); 

•  Step  2.  If  the  stack  is  empty,  go  to  Step  6;  otherwise 

•  Step  3.  Take  the  last  box  on  the  stack.  Denote  the  distance  between  the  center  of  the  box  and  the  field  point 

(i.e.,  the  center  of  the  current  cell)  as  d,  and  the  size  of  the  box  as  a.  Let  p  =  If  p  <  (3,  which  implies 

that  the  box  is  far  enough  away,  compute  the  contributions  from  all  cells  contained  in  that  box  using  only  the 
attributes  of  the  box  and  remove  the  box  from  the  stack;  otherwise 


•  Step  4.  If  the  box  contains  only  one  cell,  compute  the  contribution  of  the  cell  directly  using  the  integral  formula 
and  remove  the  box  from  the  stack;  otherwise 

•  Step  5.  Remove  the  box  from  the  stack,  put  all  its  children  boxes  on  the  stack,  and  go  to  Step  2; 

•  Step  6.  If  the  field  evaluation  at  centers  of  all  cells  is  completed,  go  to  Step  7;  otherwise  go  to  Step  1  and  start 
field  evaluation  for  the  next  cell. 

•  Step  7.  End 

Figure  1  illustrates  the  idea  for  a  2D  8  by  8  grid,  but  the  same  idea  applies  in  a  3D  grid.  To  evaluate  the  field 
at  the  center  of  cell  A,  we  need  only  calculate  the  contributions  from  the  solid  line  boxes.  Since  in  a  iV-cell  grid, 
the  field  evaluation  at  one  point  involves  O(logN)  boxes,  the  work  to  evaluate  the  field  at  centers  of  all  N  cells  is 
O(NlogN). 


Figure  1.  An  8x8  grid  of  cells  divided  into  boxes  suitable  for  evaluation  of  the  demagnetizing  field  at  the  center 
of  cell  A. 


2.4.  Multipole  Approximation  to  the  Dipole  Formula 

We  now  derive  the  attributes  for  each  box.  Assume  that  a  box  is  centered  at  the  origin  and  that  a  cell  with 
magnetization  M  is  contained  in  the  box.  Let  the  center  of  the  cell  be  ro-  From  (7),  the  demagnetizing  field  at  r, 
produced  by  the  cell,  is  approximately 


Hd(r  -  r0) 


3[M  •  (r  —  r0)](r  —  r0)—  ||r-r0  ||2M 
II  r  —  r0  ||5 


(9) 


Taking  Taylor’s  series  expansion  to  the  first  order,  we  have 


Hd(r-ro)  *  Hd(r)  +  ^^(-r0)  (10) 

3(M-r)r-M  3[fMT  +  MfT  +  (M  •  r)(I  -  5ffT)]r0 

IFF  ’  (  j 

where  I  is  the  identity  matrix.  If  we  define  the  approximation  error  in  (10)  in  a  similar  way  as  in  (8),  it  is  expected 
to  be  0(||  ro  ||2  /  ||  r  ||2).  Therefore  when  the  field  point  is  far  enough  away  from  the  center  of  a  box  (comparing 
with  the  size  of  the  box),  the  total  approximation  error  incurred  in  (5)  and  (10)  will  be  small  enough.  The  first  term 
of  (11)  is  linear  in  M,  and  the  second  term  is  linear  in  elements  of  the  matrix  Mr^.  Having  noticed  that  M  and  ro 
are  all  the  useful  information  associated  with  the  cell,  we  are  now  ready  for  deriving  the  attributes  for  each  box. 


Suppose  a  box  b  has  k  cells,  each  with  center  position  ro(s)  and  magnetization  M(s),  s  =  1,  •  •  • ,  k.  There  are  two 
attributes  associated  with  box  b:  Ai(fo)  =  M(s)  and  A 2(6)  =  X!s=i  M(s)rg  (s). 

Let  Djj .  ( i,j  =  1,2,3)  be  the  i-th  row  j-th  column  element  of  dHd(r)/dr,  and  write  Dtj  as  dfj M.  Note  that  Oij 
is  a  function  of  only  r.  Let  ©j  =  (On,  9a,  Oi 3),  *  =  1,2, 3,  and  let 

h(6)  =  [frace(©f  A2(6)),  trace(0^A2(&)),  trace(©jA2(fo))]T.  (12) 

Then  it’s  easy  to  verify  that  the  demagnetizing  held  at  r  produced  by  all  cells  inside  box  b  is  approximately 
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^Hd(r-r0(s)) 
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(13) 


The  work  to  compute  the  attributes  for  all  boxes  is  O(NlogN)  since  each  cell  is  involved  in  0{logN )  boxes.  After 
we  obtain  the  attributes,  the  held  evaluation  at  centers  of  N  cells  is  of  O(NlogN)  as  mentioned  in  Subsection  2.3. 
Thus  the  total  work  of  the  algorithm  is  still  O(NlogN). 


3.  NUMERICAL  RESULTS 

Three  methods  of  evaluating  the  demagnetizing  held  are  compared:  direct  pairwise  evaluation  using  the  integral 
formula  (Integral  Algorithm,  abbreviated  as  IA),  direct  pairwise  evaluation  using  the  dipole  formula  (Dipole  Algo¬ 
rithm,  abbreviated  as  DA),  and  the  fast  hierarchical  algorithm  (abbreviated  as  FA).  We  take  the  result  of  IA  as 
the  true  value,  and  calculate  the  root  mean  square  (RMS)  errors  for  DA  and  FA.  Let  ,  H)^t:  be  the 

demagnetizing  held  at  the  center  of  cell  ( i,j ,  k)  calculated  using  IA,  DA,  FA,  respectively.  Then  the  RMS  errors  for 
DA  and  FA  are  dehned  as 
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Table  1  compares  IA,  DA  and  FA  for  an  8  by  8  by  16  grid.  All  the  numerical  experiments  reported  in  this  paper 
were  done  on  Ultra  10  workstations  of  Sun  Microsystems. 


Table  1.  Results  of  IA,  DA  and  FA  for  an  8  by  8  by  16  grid.  (3  is  the  threshold  value  in  FA. 


Algorithm 

Time  (sec.) 

RMS  error 

IA 

89 

0 

DA 

6 

1.5  x  10-1 

FA,  (3  =  0.2 

17 

5.9  x  10~4 

FA,  [3  =  0.4 

5 

1.7  x  10"2 

FA,  [3  =  0.6 

2 

6.6  x  10~2 

FA,  [3  =  0.7 

1 

1.4  x  10”1 

From  Table  1,  we  can  control  the  accuracy  of  FA  by  varying  the  threshold  (3 .  The  algorithm  can  be  as  accurate 
as  desired.  FA  outperforms  DA  in  both  computing  time  and  RMS  error  within  a  wide  range  of  (3. 

Table  2  compares  I A  and  FA  for  grids  of  different  sizes. 

From  Table  2,  it  is  clear  that  the  computation  time  of  IA  increases  with  0(N2),  while  that  of  FA  increases 
with  O(NlogN).  When  N  goes  bigger,  we’ll  get  more  out  of  the  fast  algorithm.  Note  that  the  error  in  Table  2  is 
acceptable  in  consideration  of  the  error  incurred  by  the  discretization. 


Table  2.  Results  of  DA  and  FA  for  different  grid  sizes.  (3  =  0.4. 


Grid  size 

Time  (IA)  (sec.) 

Time  (FA)  (sec.) 

RMS  error  of  FA 

6  by  6  by  12 

15 

2 

1.3  x  10~2 

8  by  8  by  16 

89 

5 

1.5  x  10~2 

12  by  12  by  24 

928 

28 

1.7  x  10~2 

16  by  16  by  32 

5,612 

80 

1.7  x  10~2 

The  fast  algorithm  has  been  used  in  a  3D  micromagnetics  and  magnetostriction  computation  program.9  Figure 
2  shows  the  H  —  M  hysteresis  curve  of  a  3D  2x2x5  grid  computed  using  the  program,  where  H  is  the  external 
field  and  M  is  the  bulk  magnetization  of  the  ferromagnetic  body. 


Figure  2.  Hysteresis  curve  of  a  2  x  2  x  5  grid  computed  using  the  3D  micromagnetic  program. 


4.  CONCLUDING  REMARKS 

A  fast  algorithm  using  multipole  approximation  is  presented.  It’s  easy  to  implement  and  yields  computation  time  of 
0(NlogN).  By  choosing  the  appropriate  threshold  value,  we  can  make  the  algorithm  as  accurate  as  desired  while 
maintaining  acceptable  computing  speed.  It  makes  micromagnetic  computation  in  three  dimensions  feasible. 


Ixi  — 


lx  3  — 


lx  5 


APPENDIX  A.  EXPRESSIONS  FOR  lxi,  lyi ,  lzi,  i  =  1  •  •  •  8,  in  (3) 

(%  +  3  )(nz-|)  ,  _  (' nv~\){nz  +  \ ) 


( nx  +  \)\]{nx  +  \)2 +  {ny  +  \)2  +  ( nz  -  \)2 

_ (ny  +  \){nz  +  \) _ 

inx  -  \)\J{nx  -  \)2  +  (n„  +  \)2  +  (nz  +  \)2 
{ny  +  \){nz  +  |) 


?  Ijx2  — 


:5  f'XA  - 


/ - - - - :  •  —  / - ? 

(nx  +  \) \J {nx  +  \)2  +  (ny  +  i)2  +  ( nz  +  ^)2  {nx  +  \)^(nx  +  \)2  +  (ny  -  \)2  +  ( nz  -  \)2 


(nx  +  \)sj{nx  +  i)2  +  (ny  -  \)2  +  (nz  +  \)2 

_ (ny  -  \){nz  -  \) _ 

(nx  -  \)\J[nx  -  \)2  +  (ny  -  \)2  +  {nz  -  \)2 

(%  -  \){nz  -  \) 
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